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Abstract 

We implement an adaptive step size method for the Hybrid Monte 
Carlo algorithm. The adaptive step size is given by solving a sym- 
metric error equation. An integrator with such an adaptive step size 
is reversible. Although we observe appreciable variations of the step 
size, the overhead of the method exceeds its benefits. We propose an 
explanation for this phenomenon. 



1 Introduction 



Simulations including dynamical fermions remain the most challenging ones 
for lattice QCD. The standard method to simulate dynamical fermions is, 
at the moment, the Hybrid Monte Carlo (HMC) algorithm although it 
still requires large amounts of computational time. An alternative method 
to simulate the dynamical fermions is a local multiboson algorithm based on 
a polynomial approximation of the fermion matrix, proposed by Liischer0]. 
Much interest has been recently devoted to this algorithm 0] to make it as 
efficient as the HMC algorithm. 

The HMC simulation combines Molecular Dynamics (MD) evolution with 
a Metropolis test. In order to obtain the correct equilibrium, the integrator 
used in the MD evolution must satisfy two conditions; it must be: 

• time-reversible 

• area preserving 

One such integrator satisfying these conditions is the leapfrog integrator, 
which is normally used in the HMC simulation. Errors of the leapfrog inte- 
grator start with 0{At^), where At is the step size of the integrator. These 
errors cause violation of the conservation of the total energy, which must be 
corrected by the Metropolis test at the end of the MD trajectory. Let AH 
be the energy violation: 

AH = Hend — Hhegin (l) 

where Hbegin{Hend) is the total enery at the beginning (end) of the MD tra- 
jectory. The Metropolis test accepts a new configuration with a probability 
P 

prob- 

Pprob oc min(l, exp(-Ai7)). (2) 

In order to maximize acceptance of the Metropolis test, it might be preferable 
to use a higher order integrator |Q. However higher order integrators do not 
appear practical for lattice QCD since they require more arithmetic opera- 
tions (force evaluations coming from the fermionic action) than the simplest 
low-order integrator, and this overhead exceeds the gain in step-size. 
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So far conventional HMC simulations have been performed with a fixed 
step size At during the MD simulation. The local integration error does not 
remain constant in this case. When the trajectory approaches an energy bar- 
rier {Seff = —LogdetD large), it is repelled and bounces off. The curvature 
of the trajectory increases, and with it the integration error. This situation 
becomes more pronounced at small quark mass, since the height of the en- 
ergy barriers diverges in the presence of zero modes. Therefore we expect a 
behavior of the MD trajectory similar to Fig.l. Varying the step size adap- 
tively, keeping the local error constant, may be a good way to obtain a better 
integrated trajectory and may result in higher acceptance. Naively it would 
seem that this can be accomplished by calculating a local error at {p,U), 
where p = (pi,p2,---) and U = (f/i, ...) collectively represent momenta 
and link variables respectively, and then by keeping this local error constant. 
However this naive scheme is not applicable for the HMC simulation because 
it violates reversibility. 

Recently an adaptive step size method with reversible structure was pro- 
posed by Stoffer0. He constructed a symmetric error estimator which gives 
the same error value at a reflected point. The step size is then determined 
at every integration point by demanding that the symmetric error estimator 
remain constant. Stoffer implemented his method for the Kepler problem 
and obtained better results than the conventional ones. The possibility to 
apply this adaptive step size method to the HMC algorithm was stated in 
Ref . In this letter, we implement this method for the HMC simulation and 
examine its cost and its efficiency. 

2 Construction of Adaptive Step Size 

Here we construct an adaptive step size compatible with a time-reversible 
integrator. We follow the idea proposed by Stoffer 0. 
Let H be the Hamiltonian of our system, 

H=Ij:p1 + S{U). (3) 
where pi are momenta, U are gauge links, and S{U) consists of a gauge part 
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Sg{U) and a fermionic part Sf{U), 

SiU) = S,{U) + Sf{U), (4) 

S,iU)=PY.i^-j^ReTrU,ia,), (5) 

Sf{U) = <p\DD^r'<j>, (6) 



where Nc is the number of colors, is a pseudo-fermion vector and D = 
1 + kM is the Wilson fermion matrix, with k the hopping parameter. 

Call T{At) a one-step integrator. It maps momenta and link variables 
{p,U) onto {p',U'), 

T{At):{p,U)^{p',U'). (7) 
If this one-step integrator is reversible, then it satisfies 

T{~At):{p',U')^{p,U), (8) 

In this study, we use the leapfrog integrator as our one-step integrator. In 
terms of the time evolution operators]^, |^, the one-step integrator T{At) 
can be written as 

T{At) = ea;p(^L(i ^p^))ea:p(AtL(5(f/)))ea:p(^L(l 5:^,^)), (9) 

where L{») is the linear operator which is given by the Poisson bracket]^]. 
The one-step integrator requires one force evaluation represented by 
exp{AtL{S (U))) . The fermionic part of the force depends on the solution 
of a linear equation of the type Dx = (p, which is obtained iteratively at 
great expense of computer time. Thus force evaluations dominate the com- 
putation, and the cost of our algorithm can be measured in units of force 
evaluation. 

Now we define a symmetric error estimator, 

Es{p, U :At) = e(p, U : At) + e{p' , U' : -At), (10) 

where e{p, U : At) is a local error at {p, U) when the system is integrated 
by some integrator with a step size At, and the integrator maps (p, U) on 
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{p',U'). The local error is assumed to increase monotonically with At. We 
will define the local error later. If the integrator is reversible, eq.(|TDp is 
obviously symmetric under the exchange: 

{p,U,At)^{p',U',-At). (11) 

Namely, this means 

Es{p,U:At) = Es{p',U':~At). (12) 

The adaptive step size is then determined by solving a symmetric error equa- 
tion, 

Es{p, U : At) = tolerance (13) 

The tolerance should be kept constant during the MD simulation. The adap- 
tive step size determined by eq. (|1^) takes the same value at the reflected point 
(— p', U'). Therefore we flnd that an integrator with the adaptive step size 
determined by eq.(|13]) is reversible. 

Any local error can be deflned provided that it increases monotonically 
with At. Our local error is deflned as follows. 

First, we integrate [p, U) by the two-step integrator T^(At) and the one- 
step integrator T(2At). 

T\At):{p,U)~^{p',U') (14) 

T{2At):{p,U)^{p',U') (15) 

If At is not too large, (p', U') and (p', U') should be close to each other. We 
deflne the local error at {p, U) by 
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e(p, U:At)= ^(1 - ^ReTTU'l{x)U'^{x))/{W) (16) 



where V is the volume of the lattice. One could also use the momenta in the 
deflnition of the local error. Similarly, we integrate {p', U') by T^(— At) and 
T(— 2At) in the inverse time direction, 

T\-At) : {p\U') {p,U) (17) 
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T{-2At):ip',U')^{P,U). (18) 



Since the integrator is reversible, the calculation of eq.([17D is not needed. 
The local error at {p', U') is also defined like eq. ([l6D , 

e{p\ U' : -At) = £(1 - ^ReTrUl{x)U,{x))/{AV) (19) 

In the case of the leapfrog integrator of eq.(||), we need 4 force evaluations to 
construct the symmetric error estimator Es, instead of just 2 for the evolution 
T2(At). 

Eq.(|l^) is a non-linear equation. In general, it should be solved numeri- 
cally, eg. by iterative bisection. With our definition of the symmetric error 
estimator, however, we can anticipate the scaling behavior of eq.(^) and use 
it to accelerate convergence. The vector potentials evolved by the leapfrog 
integrator have O(At^) errors, 

A'^ = A'^ + 0{At'). (20) 

Therefore, 

f/'tf/' = exp{-tA'^X^)exp{tA'^X^) (21) 

^ 1 + iCaXaOiAt^) + da,pXaXpO{At^) (22) 

where are SU(3) generators, and Ca and dap are some real constants whose 
explicit values are not important here. Taking Real and Trace of eq.(|2^) 
and substituting it to eq.(0) and (|T^, we find that in the leading order 
the symmetric error estimator starts with O(At^). This behavior is verified 
numerically, as illustrated in Fig. 2: if At is not too large, the symmetric error 
estimator behaves like Eg oc At\. This property is used for solving eq. (|l3|) . 
Choose some initial At^i for the step-size and calculate Esi = Es{AtAi)- If 
Atj^i does not satisfy the symmetric error equation eq.([T3|) then input the 
second trial value At^2i which is the solution of 

tolerance .AtA2x fr>o\ 

^09{ ^ ) = Q^og{-—). (23) 
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If further trials are necessary, the following approximation can be used 0, 



At^ ^ log{At^,/At^,) tolerance 



This recurrence is continued until eq.(|T^) is satisfied to sufficient accuracy, 
and then a new configuration (p', U') is stored. Note that 2 strategies are 
available, just like for the stopping criterion of the linear solver in the force 
calculation: either the initial guess AIai is invariant under time-reversal (eg. 
it is equal to the average step-size), and the accuracy to which eq.(|^) must 
be satisfied can be set arbitrarily low; or the initial guess makes use of past 



information (eg. it is equal to the previous step-size), and eq.(|l^) must be 
satisfied exactly. We use the second method, and take the previous result as 



our initial guess. We then solve eq.(|T^) to within 5%. Since we do not solve 



eq.(|T^) exactly, we introduce a tiny, controllable source of irreversibility in 
the dynamics: the step size under time-reversal could be different by about 
5/6 ~ 1%. For this exploratory study, we have not considered this aspect 
further. 



3 Efficiency 

We examine the method for full QCD [Nc = 3) at several parameters 
{P, K, volume) hsted in Table 1. We chose /3 = in several instances, to 
eliminate the gauge part of the action and hopefully be more sensitive to the 
energy barriers coming from the fermionic part. The effect of changing the 
quark mass can be obtained by comparing cases A and B, that of changing 
the volume by comparing A and D. The adaptive step size is determined by 
the symmetric error equation eq.(|T^), with the tolerance set as per Table 1. 
The adaptive step sizes are summed up from the beginning of the trajectory 
and when the total trajectory length becomes greater than the trajectory 
length of Table 1, a Metropolis test is performed. 

Fig. 3 shows histograms of the adaptive step size at /3 = 0.0, k = 0.215 and 
0.230 ( cases A and B in Table 1 ). The distribution remains strongly peaked. 
This is also true of cases C and D. The average step size < At^ > and its 
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relative variance are summarized in Table 2, where the relative variance a is 
defined by 

As the quark mass rUq decreases, the energy barriers in phase space become 
higher and sharper, so that one would expect the variation of the step size 
to increase. Indeed this is what happens, and the relative variance in cases 
A and B increases roughly like l/m^, where rUq oc {k~^ — k~^) and Kc — 
1/4 at 13 — 0. On the other hand, as the volume increases, the relative 
variance seems to decrease sharply, like l/\/V or faster (see cases A and D). 
Perhaps this can be explained by considering the relative fluctuations of the 
effective action —Log det D: as the volume increases, the relative fluctuations 
decrease, so that the system tends to stay at some average distance from the 
energy barriers, rather than bouncing off them. 

Prom the schematic picture of Fig.l, it is expected that the approach of 
an energy barrier causes a reduction in the adaptive step size, and at the 
same time an increase in the number of iterations taken by the solver to 
converge. Fig. 4 shows At^ versus the number of iterations in the solver: the 
expected anti-correlation between them can be observed, and becomes more 
pronounced as the quark mass is reduced. 

In order to compare the adaptive method with the conventional HMC 
algorithm, we define the efficiency of the adaptive method as follows. 

First, find the fixed step-size AtHMC of the HMC simulation which gives 
the same acceptance as the adaptive step-size method. The total trajectory 
length of the HMC simulation is set to the average total trajectory length 
<traj. length> of the corresponding adaptive step size method's case. We 
performed the HMC simulations with several step-sizes and determined the 
corresponding step-size of the HMC algorithm by interpolating those results. 
The results of the corresponding step-size are summarized in Table 2. For 
the acceptance of the adaptive step-size method, see Table 2. 

Then, define the gain by 

QA =< AtA > /AtHMC- (26) 
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When gA > I, the adaptive step-size method really takes larger steps on 
average, without compromising the acceptance. However the real efficiency 
of the method can only be assessed by taking into account the overhead 
of determining the adaptive step-size, since additional force evaluations are 
necessary. 

From the definition of our symmetric error estimator eq. ([I^ -(p!9D, we 
know that one construction of Es needs 4 force evaluations. Call Rt the 
average number of trial steps needed to solve eg. (p!3D . After ARt force evalu- 
ations, we have determined the step-size At a and use the integrator T^(At^) 
(eq.|T^) to advance the dynamics by 2 steps At a- Therefore the cost of the 
method is 2Rt force evaluations per step, compared with 1 force evaluation 
per step for standard HMC. 

Real efficiency will be achieved if the number of force evaluations per unit 
time decreases, ie. if 2Rt < Qa- Results for the gain qa and the cost 2Rt 
are summarized in Table 3. For all cases we studied, real efficiency is not 
achieved. 

It is disappointing to see how small the gains qa are. The reason for such 
small gains can be understood by considering the behavior of the Hamilto- 
nian. The dependence of the energy violation at each step with the step-size 
is, in general, non-linear. Therefore it is not necessary that a small local er- 
ror correspond to a proportionally small energy violation of the Hamiltonian. 
Fig. 5 shows |AiJ|, the absolute value of the energy violation after one inte- 
gration step, versus the local error Es- The 2 clusters of points correspond 
to fixed step sizes At = 0.04 and 0.08. No strong correlation between Es and 
AH can be observed. This is further evidenced by the dashed lines, which 
are the result of fitting to a scaling law v<AfP> = E^^. for the larger 
step-size, h is almost zero. Therefore, it becomes clear that fixing Es and 
varying At adaptively cannot have a strong effect on the acceptance, which 
solely depends on AH . 

Two approaches could be used to improve the efficiency of our scheme: 
i) decrease the overhead: instead of estimating the error by comparing T^(At) 
with T(2At), one could replace the latter by an Euler integrator, which 
requires no additional force evaluation. Note that the error eq.(|l^) remains 
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symmetric under the exchange {p, U, At) < — > {p', U', — At) even though the 
Euler integrator is not time-reversible. The problem we found with that 
approach is that, for the large step-sizes used on our small lattices, the error 
eq. ([TO|) no longer obeyed a simple scaling law eq.(p2t) as a function of the 
step-size. Then the number of iterations needed to solve eq.(|l3D increased, 
defeating the expected reduction in overhead. On larger lattices with smaller 
step-sizes, this problem would be milder. 

ii) change the definition of the error eq. ([T9|) , so that it is better correlated 
with |Aif|, the energy violation at each step. Note that \AH\ itself cannot 
be chosen, because it does not increase monotonically with the step-size: in 
that case eq.(0) admits multiple solutions; the overhead of converging to one 
of them, and the same one under time-reversal, increases considerably. With 
our definition eq . ([19|) , the error is only weakly correlated with |AiJ|, but the 
situation again seems to improve with smaller step-sizes, on larger lattices 
(compare the 2 dashed lines in Fig.5). Nonetheless it would be desirable to 
control the step-size with a more relevant quantity than eq. ([T9|) , since all 
that matters in the end is energy conservation. 

Finally, instead of varying the step-size, one could vary adaptively the 
couplings of the Hamiltonian H, eq.(0), at each step, or even include some 
new operators in H, trying to tune them so as to best conserve energy. The 
general difficulty with that approach is to find an error eq. (|l9|) which varies 
monotonically with the couplings of H. 

4 Conclusion 

We have implemented an adaptive step-size method for Hybrid Monte Carlo 
simulations, and tested it at several parameters [jS, k, volume). The relative 
variance of the step-size increases for small quark masses and small volumes. 
The average step-size seems somewhat larger than the corresponding fixed 
step-size at same acceptance. But this gain is more than offset by the over- 
head of determining the adaptive step-size. It seems very difficult to achieve 
real gains in efficiency, because conservation of energy, which is necessary for 
high Metropolis acceptance in the HMC algorithm, is poorly correlated with 
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the conventional error governing the adaptive step-size. 

A plausible extrapolation from our results would indicate that the relative 
variance of the step-size scales like m~^V~^^'^, ie. as (niT^L)''^, where is 
the pion mass and L the physical size of the lattice. This quantity normally 
remains constant as the continuum limit a — > of the lattice theory is taken, 
so that the relative fluctuations in the adaptive step-size At would tend to 
a constant. Even if this analysis is no more than plausible at this stage, it 
is clear that the two limits ruq ^ 0, V ^ oo have opposite effects on the 
fluctuations of At, making it unlikely that such fluctuations become very 
large on present lattice sizes. This observation is consistent with the limited 
fluctuations (a factor 2 or so) in the number of iterations needed by the solver 
to compute the force in the largest QCD simulations 

Thus it appears that QCD is much "easier" to simulate than the Kepler 
problem: in lattice QCD, the force on the gauge links varies little in magni- 
tude, and the curvature of the Molecular Dynamics trajectory is rather small. 
One intuitive explanation is that the QCD force is dominated by short-range 
UV contributions, which drown the IR component sensitive to the energy 
barrier detD ~ 0. 

We thank D. Stoffer for helpful discussions. T.T. is supported in part by 
the Japan Society for the Promotion of Science. Ph. de F. thanks Hiroshima 
Univ. and Tsukuba Univ., especially Profs. O. Miyamura and Y. Iwasaki, 
for hospitality during this project. 
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Figure Captions 

Figure 1: Schematic behavior of a Molecular Dynamics trajectory in con- 
figuration space. detD is the determinant of the Dirac matrix. 

Figure 2: Adaptive step size At a versus symmetric error Eg, for two 
configurations of size 4^ at k = 0.230. The straight line Es oc At\ is shown 
for comparison. 

Figure 3: Histograms of adaptive step size At a, for 4^ lattices at k = 0.215 
and 0.230. The logarithmic scale shows the increase with k of the relative 
variance. 

Figure 4: Adaptive step size versus number of iterations in the solver 
(BiCG75). 

Figure 5: \AH\ versus the local error £'5, on a 4^ lattice at k, — 0.215. 
The step size is fixed at At — 0.04 and 0.08. AH is the change in the total 
ft 01 one integration step. The dotted lines result from fitting to the 
form V < AH^ > — Eg, and show the correlation (or absence of) between 
the 2 quantities. < AH^ > is obtained by dividing the data in 10 bins and 
averaging the values AH^ in each bin. 
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case 


/? size 


K 


traj. length 


tolerance (±5%) 






A 


0.0 4^ 


0.215 


0.8 


10-^ 






B 


0.0 4^ 


0.230 


0.4 


8 X 10-6 






C 


5.4 4^ 


0.162 


1.0 


10-6 






D 


0.0 8^ 


0.215 


0.3 


10-' 












Table 1: Run parameters 






case 


< At A > 


a{%) 


<traj. length> acceptance(%) 


AtHMC 


A 


0.0911(3) 


3.3 


0.91 


36(2) 




0.0897(13) 


B 


0.0431(1) 


5.6 


0.44 


57(2) 




0.0419(09) 


C 


0.0673(1) 


0.8 


1.08 


87(2) 




0.0688(80) 


D 


0.03281(2) 


0.6 


0.33 


63(2) 




0.0328(10) 



Table 2: Results of the adaptive step-size method and fixed step-size AtHMC 
of the HMC algorithm. <traj. length> stands for the average total trajectory 
length. The fixed step-size AinMC is determined so that it gives the same 
acceptance as the adaptive step size method. 
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